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Abstract. For general anisotropic linear elastic solids with smooth bound- 
aries, Rayleigh-type surface waves are studied. Using spectral factorizations of 
matrix polynomials, a self-contained exposition of the case of a homogeneous 
half-space is given first. The main result is about inhomogeneous anisotropic 
bodies with curved surfaces. The existence of subsonic free surface waves is 
shown by giving ray series asymptotic expansions, including formulas for the 
transport equation. 



1. Introduction 



Ravleighl (|1887I ) discovered waves which propagate along a traction-free plane 
surface of an isotropic and homogeneous elastic solid. The surface wave speeds 
are subsonic, that is, they are strictly less than the wave speeds of interior body 
waves. F urthermore, t he amplitudes attenuate exponentially with distance to the 
surface. Svngel ( 1957 ) r aised doubts about existence of Rayleigh-type waves in 
anisotropic solids. IStrohl (|l962l ) pointed out that these doubts were unfounded, and 
he introduced a sextic eigenvalue problem which became useful in the theory of free 
surface waves in anisotropic solids. In the early 1970's, the existence and uniqueness 
problem of free surface homogeneous plane waves in a semi-infinite half-space was 
settled by Barnett, Lothe, and coworkers. For any given horizontal propagation 
direction they showed that there is at mos t one free surface wave s peed, and they 
gave criteria for the existence of such waves. Lothe fc Barnett ( 1976[ ) rederived their 
results by the surface impedance tensor method. The surface impedance tensor 
relates the surface displa cement to the surface traction r equired to sustain it. This 
tensor was introduced bv llngebrigtsen fc Tonniii3 ( 19691) . Detailed p resentations of 
the existence and uniqu eness results were given by Chadwick fc Smith (197 7) and 
Barnett fc Lothd (|l985h . Much later, iMielke fc Ful (|2004l ) simplified some proofs of 
the Barnett-Lothe theory by using a Ricatti equation satisfied by the impedance 
tensor. A crucial property of the tensor, the positive definiteness of its real part, 
follows from an integral identity which, in the original treatments, arises somewhat 
magically by averaging over rotations in the plane spanned by the normal to the 
surface a nd the propagation dir e ction. Existence of subsonic free surface waves was 
shown by iKamotskii fc KiselevI (|2009l ) with a completely different approach based 
on the variational principle. 

Concerni ng Rayleigh-type waves in inhomogeneous elastic solids with curved 
boundaries, IPetrowskvl (|l945l ) exhibited the following locality principle: If a sur- 
face wave exists, the velocity of its discontinuity at a given point must be equal to 
the velocity in the homogeneous elastic half-space which is obtained by freezing the 
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elastic parameters at that point. This locahty principle is efficiently implemented 
by asymptotic ray methods, which substitute, in the high-frequency regime, stan- 
dard plane waves by geometrico-optical 'plane waves'. For inhomogeneous, isotropic 
elastic solids with curved boundaries, ray methods involving sums o f complex plan e 
waves were successfully app lied to t he free surface wave problem bv lBabichl(|l96ll) . 
Babich fc Rusakoval (|l962[ ). and by iKaral. Jr. fc Kelleij (|l964l ). The ray method 



works, for smooth data, under assumptions on the geometry and on the elastic pa- 
rameters which are less restrictive than those needed for finding analytic solutions. 
More importantly, salient features of high-frequency waves such as wave fronts and 
amplitudes are captured directly by ray methods. Usin g the existence and reg ular- 
ity theory of linear hyperbolic equations, as was done bv lCourant fc LaxI (|l956h . one 
can correct asymptotic solutions, withou t changing the important high-frequency 
properties, into genuine exact solutions. iGregorvi (|l97ll ) compared some analytic 
representations of surface waves and corresponding ray approximations. Rayleigh 
surface waves in inhomogeneous, anisotropic elastic bodies were studied with the 
complex ray method by Nomofiloy (1979)- The amplitude of a geometrico-optical 
wave is governed by a transport equation which, in the case of free surface waves, is 
quite complicat ed. Efforts to solve the transport equations culminated in the work 
of B abich fc Kir pichnikova (2004)) , where detailed formulas for the amplitude and 
the Berry phase of a Rayleigh surface wave were obtained. 

In the present paper we study Rayleigh-type surface waves in anisotropic elastic 
solids with smooth surfaces and inhomogeneities. We give a self-contained presen- 
tation of the theory of Barnett and Lothe, and we incorporate it into a ray theory. 
Under the same conditions as for homogeneous half-spaces, the existence of sub- 
sonic free surface waves is proved. Transport equations for leading amplitudes are 
established in a way which enables their numerical solution. 

Our approach uses spectral factorizations of the acoustic tensor which is regarded 
as a second order matrix polynomial in the variable conormal to the boundary. This 
allows to conveniently lump together the relevant complex eigenvalues and to avoid 
some cumbersome and unnecessary considerations. Moreover, the factorizations 
reduce the elastodynamic system, in the subsonic regime, to a first order system. 
In general, when the surface is curved and the solid inhomogeneous, the first or- 
der system is not differential but pseudo-differential. The zero traction boundary 
problem is transformed into a pseudo-differential wave equation on the space-time 
boundary. The associated principal symbol is the surface impedance tensor. The 
pseudo-differential wave equation is only defined microlocally over the subsonic re- 
gion. Still it can be treated ray theor e tically since it has the structure of a so-called 
real principal type system (jPencken . 119821 ). In the microlo cal ana l ysis li terature 
such an approach was carried out for the isotropic case bv iTavloil ( 1979f ). where 
the subsonic region is called the elliptic region because the theory of elliptic bound- 
ary value problems applies in the reduction to the boundary. The idea to split 
a surface wave problem into an elliptic boundary problem and a wave propaga- 
tion problem along the space-time boundary h a s also been used in related work 
( Kaplunov et al. . 2006 ) in elasticity. Nakamural (jl991 ) treated the Barnett-Lothe 
method from the p oint of vie w of microlocal analysis. The present paper uses ideas 
developed in (Han sen fc Roh rig. 2004; Hansen, 201lll201^, and it is an attempt to 
communicate the ideas and results to researchers of elastic surface waves who may 
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not be experts in microlocal analysis and its main machine, the pseudo-differential 
calculus. 

The paper is organized as follows. Section [5] recalls the free surface traction 
problem in differential geometric tensor notation. In sections [3] and |4] the subsonic 
free surface wave theory for homogeneous elastic half-spaces is redeveloped using 
division theory of matrix polynomials as the basic tool. In the appendix of the 
paper, we give a complete and self-contained presentation of those parts of division 
theory which we employ. The remaining sections deal with inhomogeneous solids 
with curved boundaries. For the benefit of readers who are not familiar with pseudo- 
differential operators, Section [5] contains an account of core properties of pseudo- 
differential calculus and its relation to asymptotic expansions. A ray theory for 
pseudo-differential systems is developed in Section |6l In Section [7] the displacement 
boundary problem is solved for surface displacements concentrated in the subsonic 
region. A factorization of the elastodynamic operator is constructed and used in 
that section. Theorem 18.11 in Section |8] is our main result on the existence of 
Rayleigh-type waves. The transport equation for the leading amplitude of subsonic 
free surface waves is treated in some detail in Section [9l 



2. Equations of linear elastodynamics 

Let i? be a linearly elastic body with density p > and stiffness tensor C = 
\^(jvkij_ -^^ijjg elasticities satisfy the standard symmetries, 

and the strong convexity property, 

(2) C*^'"eye^>0 if e,, = e.,, e = [e,,] ^ 0. 

(We use the summation convention, and we denote a point, if at all, by its coordi- 
nates X = {x^,x'^,x^). A bar denotes complex conjugation.) Since the elasticities 
are real, it suffices to assume that holds for symmetric tensors e which are 
real. Assumptions ^ and ^ say that C defines an inner product on symmetric 
2-tensors. 

The strain tensor field e measures the deformation of the geometry of B caused 
by an (infinitesimal) displacement u. The geometry is given by the Riemannian 
metric tensor G — [Gij]. Precisely, strain is the symmetrized covariant derivative 
of displacement: 

£ke = {uk-e + w<!;fe)/2, Uk-e = Ukj - ^igUj. 

Here F^^ = G'^^{Gij,k + Gkjj ~ Gke.j)/^ are the standard Christoffel symbols of 
differential geometry. We precede a coordinate index by a comma or a semicolon to 
denote a partial or a covariant derivative, respectively. The stress field, a = [cr*-'], 
represents forces. It is related to strain via Hooke's Law: 

(3) a'' - C^'^'eki = C'^^'uk-^i. 

The inverse of G is G^^ = [G"-']. By raising and lowering indices one switches 
between covariant and contravariant components, e.g., = G^^Uk and Uk — GkjU^ ■ 
The metric G defines the length element ds by ds^ — Gij dx' dx^ , and the volume 
element d^ = ^/g dx, g = det G. 
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Let S denote the (smooth) boundary surface of B, dS its surface element, and v 
the interior unit normal field. Assume that B is source-free, and that S is traction- 
free. Then the displacement field satisfies the free surface boundary problem: 

(4) pi;^-a'\^=0, a'^i^,\s=0, 

i = 1,2, 3. The differential equations are the 3x3 system of elastodynamics. The 
equations (jU arise as the Euler equations of the Lagrangian which is the functional 
given by 



C{u) = J J {pu'u, - a'^e^j) dV dt. 
The divergence of the stress tensor is 

(5) a'{^^ = a'\^ + n^a'^^ + T^y^ = g-'^' {g'^'a^') + Tlj'^"' ■ 

The last equality follows because (logg),j — ^T^j. 

The differential geometric formulation of elastodynamics, which we just recalled, 
does not depend on the choice of a particular coordinate system. Elasticities, 
strains, stresses, and displacements transform as tensors under changes of coor- 
dinates. No generality is lost when we assume that the surface S agrees, near some 
chosen point, with a coordinate plane, e.g., — 0. In the following, we always use 
coordinates which are adapted to S in the following sense: The normal coordinate 

is the signed distance to 5* (such that a;'^ < in the exterior) , and the horizontal 
coordinates x^ and x'^ are constant along the geodesies which intersect the level 
surfaces of x^ orthogonally. The horizontal coordinates are completely determined 
by their restrictions to S. The metric tensor satisfies G^^ = 1 and G^^ = G^^ = 
if A < 3. The surface area element equals dS = ^ dx^ dx^, and r* = a^^ is the 
i-th component of the traction r at a given level surface x^ =constant. 

In the study of surface waves, it is useful to separate, in the elastodynamic 
equations and in the formula for the traction, differentiations in normal direction 
from differentiations in the horizontal (and time) directions. We write 

(6) a''^^ = C*3'=3„fc,33 + (C^*'"^ + C'='^^)^ifc,A3 + C'^'^^'uu^^^, + B'^^'uk^i + B'^uk, 

(7) = C'^'^^uk^s + C'3'=^Wfc,A - C'^'^'r^uk, 

where the implied summations are restricted to A,/x < 3. The coefficients of the 
first order derivatives are 

D — O .j~ ^ ^ ran^^ ^ jm + "-^ ^ raj ■ 

We have no need to know the i?*'^'s explicitly, only their vanishing if the metric 
tensor G is constant. 



3. Surface waves in homogeneous half-space 

In this section we assume that the elastic solid is homogeneous and occupies a 
half-space in Euclidean space. Thus the density, the elasticities, and the metric 
tensor are constant. In particular, the Christoffel symbols are zero, and covariant 
derivatives are ordinary derivatives. 

We use adapted coordinates. Then x'^ > corresponds the half-space filled by 
the elastic body. The unit normal field v at S, which points into the interior, has 
the components v\ ~ = ^ \. 
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The equations of the elastodynamic free boundary problem are as follows: 

(8) pu, - Q^'^S^,, 33 _ (Q3feA ^ C'^S^A)^^ _ c/'=''ufe,AM = 0, 

(9) C,3'^'3^fc,3 + C,;'^V.A = atx^^O. 

Again summation is restricted to A, /i < 3. We have lowered the index i in order 
that the elastodynamic and the traction operators map from covariant to covariant 
components. 

Interior and reflected plane waves are found with the time-harmonic ansatz 

u{x,t) = exp(iK(^x — ct))U, = (,jX^ ^ i = \/— T, 

where k > is the wave number, c > the wave speed, ^ = [^j] the unit prop- 
agation direction, and ^/c the slowness vector, \^\ — 1. The ansatz satisfies the 
elastodynamic system if and only if the amplitude vector U — \Uk\ lies in the null 
space of the acoustic tensor [C/'^'^^j^^ — c^p5^'^]. To obtain waves which propagate 
along the surface, this ansatz is modified by allowing ^'^ to be complex with positive 
imaginary part. 

Let rj — [rjj] be a real horizontal unit vector, which means that I77I = 1 and 77 is 
orthogonal to v, i.e., 773 — 0. Consider the acoustic tensor at 77 + su: 

A{s) = s^Aq + s(Ai + Al) + ^2 - c^pl- 

Here / = [5j^] denotes the 3x3 unit matrix, and 

(10) A,^[C,^^% A,^[C,^^'m]. A,^[C/''r^,m] 

are real 3x3 matrices; Af is the transpose of Ai with respect to the inner product 
defined by G. Because of the symmetries of the stiffness tensor C, the matrices 
Aq, A2, and, for real s, A{s) are symmetric with respect to G. The wave speed 
c > is said to be subsonic if A{s) is non-singular for every real s. By the strong 
convexity of C, ^0 is positive definite, and there exists a positive limiting wave 
speed Coo = Coairj) such that c is subsonic if and only if c < Coo- Notice that A{s) 
is positive definite for real s if c is subsonic. 

We use the following time-harmonic ansatz to find surface waves which have 
subsonic horizontal slowness 77/c: 

(11) u{x, t) — exp{iK{rix — ct))U{KX^). 

The amplitude U{z) at z = kx^ > is to be determined. We introduce the 
differential operator DU{z) = —i-j^U{z). The ansatz (fTTj) satisfies equation (|5]) if 
and only if 

(12) AoD^U + {Ai + Al)DU + (A2 - c'pI)U = 
holds. 

The matrix polynomial A(s) satisfies the assumptions of Proposition I A. 1 1 in the 
appendix. Therefore, there is a unique complex 3x3 matrix Q with spectrum in 
the complex upper half-plane and such that, with Q* denoting the adjoint of Q, 

(13) A{s) ^ {si - Q*)Ao{sI ^ Q) 
holds for complex s. Now (IT2|) can be rewritten as 

(14) {DI -Q*)Aq{DI -Q)U = 0. 

An analogous formula holds with Q replaced by a matrix P, which has its spectrum 
in the lower complex half-plane. The solutions of DU ~ PU and DU = QU 
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span the space of solutions of ([T2|) . The solutions of ([T2|) . which stay bounded as 
Kx'^ — > +00, are precisely the solutions of DU = QU, which are given by U{z) = 
exp{izQ)U{0). Moreover, as z +00, these solutions decay exponentially. 
The surface traction of the waves pT|) just constructed equals 

T = iK{AQDU{0) + AiU{0)) = -kZU{0). 

Here Z is the 3x3 surface impedance tensor, 

(15) Z=-i{AoQ + Ai), 

first introduced by llngebrigtsen fc Tonningj (|l969t) . Summarizins;, we have found 



that, precisely when ZU{0) = holds, the ansatz ([TT|) leads to a solution of (|8]) 
and ([9]), which decays into the interior. Since Q depends smoothly on the subsonic 
wave speed c and on 77, so does Z . Suppressing the dependence on 77, we often write 
Z(c) for the impedance tensor. 

In the early 1970's the problem of uniqueness and existence of subsonic surface 
waves in homogeneous half-space was completely solved. The following theorem, 
which we reprove in the next section, states the solution in terms of the impedance 
tensor. 

Theorem 3.1 (Barnett & Lothe). For a given propagation direction, there exists 
at most one subsonic free surface wave. A subsonic surface wave exists, if and only 
if det Z{c) < holds for some c in the range < c < Coo- The wave speed is the 
unique solution, in the range < c < Coo, of the secular equation det Z(c) — 0. The 
surface displacement of a subsonic surface wave belongs to the null-space of Z{cfi). 
The null-space of Z{cfi) is one- dimensional and contains no non-zero real vector. 

4. The surface impedance tensor 

The surface impe dance tensor method of subsonic surface wave theory, devel- 
oped bv llngebrigtsen fc Tonnindl ( 19691) and lfiarnett fc Lothe ( 1985( ). relies on the 



following properties of Z. 

Proposition 4.1. For < c < Coo the following hold: 

(a) Z{c) is Hermitian. 

(b) Zifi) is positive definite. 

(c) The real part of Z{c) is positive definite. 

(d) At most one eigenvalue of Z{c) is non-positive. 

(e) The derivative dZ{c)/dc is negative definite. 

(f) The limit Z{coo) = linicfcoo ■^^(c) exists. 

Theorem 13.11 follows from Proposition 14.11 combined with arguments of the pre- 
vious section. Indeed, it follows from the proposition that the derivative of the 
determinant with respect to c is negative at its zeros. Hence det Z{c) = has at 
most one zero. Since detZ(0) > 0, a zero c = cr exists if and only if the determi- 
nant becomes negative for some subsonic c. The assertions about the null-space of 
Z{cr) follow from (jcj) and 0. 

No tice that det Z{c) < holds if and only ii cr < c < Coo- Barnett fc Lothd 
()l985 . Theorem 12) state an existence criterion in terms of Z(coo). 



Proof of Proposition \4.1\ Following iMielke fc Ful (|2004l ). we use the Ricatti equa- 
tion 

(16) {Z-iAj)AQ\Z + iAi)=A2-c^pL 
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Observing that Q = iv4g ^(Z + iAi), ([T6|) is seen to be equivalent to 
Aog' + (Ai + Al)Q + A2- c^pl - 0, 

which is the solvency equation ((50)) of the factorization (IT3l) . 

Passing to adjoints, one recognizes that (fT6|) also holds when Z is replaced by 
Z*. Subtracting the two Ricatti equations, one obtains 

Q*{Z-Z*)-{Z-Z*)Q = 0. 

Because Q and Q* have disjoint spectra, this Sylvester equation is non-singular. 
Hence Z* = Z, proving (jsj). 

To prove (gj), differentiate with respect to c, and get 

iQ*Z -iZQ = -2cpl, Z=dZ/dc. 

This Lyapunov-Sylvester equation has the unique solution 

/■oo 

Z = —2cp / exp{—isQ*) exp{isQ) ds, 
Jo 

which is negative definite. 

Next we prove the positive definiteness of Z{0) ~ [Z^^]. Assume that w is 
a complex vector such that Z^^WkWl < 0. We have to show that w — 0. The 
exponentially decaying solution of (|T2|) with surface displacement e"'^w is 

u = e*''^ exp{ix^Q)w. 

Denote by e = [cke], (-u = {diUk + dkUi)/2, the strain tensor and by r = —Z{Q)w 
the surface traction of the (complex) displacement field u. Integrate the divergence 

over the half-line > 0, and get the energy identity 

J C'''^e^Je^^dx^ ^ -T^w^. 

By assumption, the right-hand side is non-positive. Hence, in view of the strain 
vanishes on the half-line > 0. In particular, 

en = irjiui, £12 = i(?7iU2 + ?72Ui)/2, and £33 = iiQu)^ 

all vanish. Without loss of generality, we assume that 771 7^ 0. It follows that 
m — U2 — 0. Moreover, i93|u3p = 2Re {i{Qu)3U3) = 0. Since U3 tends to zero as 
x"^ — > 00, this implies U3 = 0. Therefore w = 0, which proves assertion (jb|) . 

Denote by A{s) the polynomial The integral formula for Q, of the 

appendix, implies 

iZ I A{s)-^ds= I {sAn+ Ai)A{s)-'^ ds, 

if the closed contour ^b. consists, with i? > sufficiently large, of the interval 
[-R,R] and the arc i?e"^, Q <Lp <'k. Notice that AnA{s)-^ = s^^/ + 0{s-^) as 
|s| — > 00. Therefore, letting i? — >■ cxd, we obtain 

/oo /"OO 
A{s)~^ ds ^ nil + {sAo + Ai)A{s)-^ ds, 
-00 J —00 
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the integral on the right being a principal value integral. The integrals are real 
matrices. Moreover, the integral on the left-hand side of ([T7|) is positive definite. 
This proves the assertion 

If (jd| were not true, then, for some velocity c, the impedance tensor Z{c) — [Z^^] 
would have only one positive eigenvalue. Furthermore, its eigenspace would be one- 
dimensional. Since space is three-dimensional, we could then choose a real vector 
w orthogonal to this eigenspace. But then Z^^WiWj < 0, contradicting the 
positive definiteness of the real part of Z{c). Hence holds. 

To prove U we first derive a bound on the operator norm ||^(c)|| of Z(c). Notice 
that ||.^(c)|| equals the maximum of the absolute values of the eigenvalues of Z{c). 
By (jej), Z{c) < Z{0) holds with respect to the usual ordering of Hermitian matrices 
by the cone of positive definite matrices. Therefore the eigenvalues of Z(c) are not 
greater than ||Z(0)||. The sum of the eigenvalues of Z{c) is positive because it is 
equal to the trace of the positive definite matrix Re Z{c). Therefore the modulus of 
a negative eigenvalue, if it exists, is less than the sum of two positive eigenvalues. 
Hence we have ||-^(c)|| < 2||Z(0)||. By compactness, there exist limit points of 
Z{c) as c t Coo- Using (je| again, we find that Z^, < Z{c). It follows that all limit 
points are the equal, which implies assertion □ 

Remark 4.1. Formula (|17p is a variant of the Barnett-Lothe integral representation. 
In fact, passing to the adjoint of (|17|) and making a substitution (p — cot(s) to 
replace the integration variable s by an angle (/>, one rederives the formula (2.19) of 
Barnett fc Lothel (|l985l ). 



Remark 4.2. If the coefficients of the matrix polynomial depend smoothly on pa- 
rameters, then also the wave speed cr depends smoothly on these parameters. Since 
the zeros of the determinant are simple, this follows from the implicit function the- 
orem. 



5. Asymptotic expansions and pseudo-differential operators 

If the elastic body is inhomogeneous and the boundary surface curved, we do 
not expect that an exact operator factorization ([14]) holds. However, using pseudo- 
differential operators, we shall factorize the elastodynamic operator up to negligible 
errors, and we shall use this to exhibit asymptotic subsonic surface waves. This 
will be detailed in the following sections. As a preparation, we summarize stan- 
da rd results on asymp t otic expansi o ns an d on pseud o-differential opera t ors. Refer 



to iHormandeij (|1965D . iHormandeil (|l985l . 18.1), or lAhnhac fc GerardI (|2007f) for 



expositions of the basic pseudo-differential calculus. 

We consider vector- valued functions u{x] uj) of points x in d-dimensional space 
which oscillate rapidly as the frequency parameter uj tends to -\-oo. More precisely. 



(18) 



u{x; uj) 



k=0 



where the series is an asymptotic series in the space of smooth functions x. (Often 
we simply write = instead of ~ , although the equation may be true only in the sense 
of asymptotic expansions.) The phase function 9(x) is real-valued with derivative 
0'{x) ^ 0. We use the operators Dj = —idj, where dj = d/dx^ denotes partial 
derivative with respect to the j-th coordinate, x^ . By the Leibniz' product rule. 
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/ = DjU has the asymptotic expansion 

(19) fix; w) ~ c."e-^(-) Uu:)-^F^k{x) 

with m = 1 and top order coefficient Fq{x) = dj9{x)Uo{x). A linear differential 
operator order m, A, is a polynomial in Dj of order m with smooth coefficients. 
Also / = Au has an expansion ([T9|. 

The pseudo-differential calculus assigns to certain functions A{x,^), which are 
defined on 2(i-dimensional phase space, operators A = A{x, D) by 

(20) (Mix) = (27r)-'' / / e'(^"^)-«A(x,<e)u(y) dy d^ 



Here x and ^ denote position and (generalized) momentum variables, respectively. 
By the Fourier inversion formula, A is the identity when A = 1. The function 
A{x, 5) is called the (full) symbol of the operator A. We assume that the symbol A 
belongs to a standard class — of symbols of order at most to. Moreover, the 
symbol A admits, as |^| — )■ oo, an asymptotic expansion in homogeneous functions: 

A{x,0 - E.to'^-^'^'^'^^' ^-^■(^'*^) = t™-^ A_,(x,0 for t > 0. 

(Strictly speaking, homogeneous functions are, unless they are polynomials, not 
symbols. This technicality is overcome by multiplying with cutoff functions which 
insure smoothness of symbols at ^ = 0.) The symbol A{x,£^) and the associated 
operator A are said to be of order m if the principal symbol Aq{x, ^) is not identically 
zero. Symbols may be matrix- valued. Rather than using A^i{x,^) it is preferable 
to work with the subprincipal symbol Asuhix,^), 

(21) A.ub = A-i-ll d^Ao/dx'd^j. 

A pseudo-differential operator is a differential operator if and only if its symbol is 
a polynomial in the ^ variable. In particular, the symbol of Dj is ^j. 

Every pseudo-differential operator of order m has the property that it maps 
an asymptotic sum (fTSi) into an asymptotic sum (1191). Moreover, th is property 



is characteristic of pseudo-differential operators; see iHormander (Il965l) . bmce we 



need explicit expressions for Fq and F-i we indicate how (|T9|) is proved using the 
method of stationary phase. The method is applied to the summands of (fT8|) : 



oi-^)(AUe''''^)ix) = {u:/2nf / / e*"("-^)-«-*"(^(")-^(^»A(x, wOC^(y) dy d^ 



The stationary point of the phase 

<i> = (x-y).c-(0(x)-%)), $; = o = $^, 

is at y = ^ = 0'{x). The inverse of the Hessian matrix of <f> at the stationary 
point is the 2c? x 2(i-matrix 

"0 / 

/ e"{x) 



where H{x) = 



On a formal level, the stationary phase expansion is given by 

e-'-^'>(^){AUe'^'>){x) - exp (^{iij)-\dy ■ + h"dy ■ dy)^ (e^^^C/), 
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where the exponential of differential operator s has to be replac ed by its formal 
Taylor series to give the asymptotic expansion (jHormandeii Il990l Theorem 7.7.5). 
Here p is the remainder of the second order Taylor expansion of the phase, it 
vanishes to third order at the stationary point where the expressions have to be 
evaluated. Thus the principal coefficients in ((T9)) are 

Fo{x)^Ao{x,e'{x))Uo{x), 

and, suppressing the arguments x and ^ — 9'{x) in writing, 

P-^-Y.^^ + -,[Y.o^.9^^)u^ + ^^-^Uo + AoU^,. 
Observe that the double sum enclosed in parenthesis equals 

Therefore, 

For k — 2,3,... there are formulas for F^^, which differ from the formula for F^i 
by a shifted index and by an additional term which is a sum of (derivatives of) 
U-j with j < k — 1. The additional term arises from the higher order terms in the 
stationary phase expansion. We emphasize that the asymptotic expansions and the 
formulas for the coefficients F^k hold for a general system A of pseudo-differential 
operators. 

The composition C — BA of pseudo-differential operators A and B is again a 
pseudo-differential operator. There is a formula for the asymptotic expansion of the 
full symbol C{x,S,). On the principal symbol level, symbols multiply: Co = BqAq. 
On the next level, 

(23) C_i = Bo^-i + S-i^o -i^X5c,^o)(9.Mo) 
holds. Define the Poisson bracket 

{P, Q} ^ Y. -^^Pl^^o){dQ/dx^) - {dP/dx^){dQldij), 

where the symbols need not be scalar but can take square matrices (of equal di- 
mensions) as their values. It follows that the composition of operators is given on 
the principal and on the subprincipal level by 

(24) Co = BqAo, Csub = So^sub + SsubAo + ^{Bo, Ao}. 

The (formal) adjoint of A with respect to a given scalar product is also a pseudo- 
differential operator, _B, and its principal and subprincipal symbol are given as 
follows: 

Bq = Aq, Bsuh = ^sub- 

The star denotes the (Hermitian) adjoint. (The symbol formulas arise more nat- 
urally if Weyl quantization is used instead of the Kohn-Nirenberg quantization, 
j4 A as in (f20| . which we are using here.) 

Under ellipticity assumptions, one constructs inverses, square roots, and powers 
of pseudo-differential operators which are again pseudo-differential operators. The 
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error or remainder terms in such constructions are typically negligible operators 
with symbols belonging to S^°° = DmS™. In particular, two pseudo-differential 
operators with symbols, which are equal when |^| > 1, differ only by a negligible 
operator. Negligible operators map distributions into smooth functions. If u has 
an symptotic expansion (|18p and if A is a pseudo-differential operator which is neg- 
ligible in conic neighbourhood of the set of {x,uj6' (x)), uj > 0, then the asymptotic 
expansion oi f — Au is trivial, i.e., all coefficients F-k — 0. 

The operator L of elastodynamics is a second order differential 3x3 system, 

(25) {Lu)t = pili ~ (7-^.^. 
The principal symbol Lq of L is the acoustic tensor, 

(26) Lo{x,0 = [C/''H^)^j^i-peo5,H^)]- 

Here we included the time coordinate as = t, and the dual variable (frequency) 
as ^0- The implied summation extends over the spatial indices, i.e., j,£ > 1. 
The pseudo-differential operator which will be most important to us is the surface 
impedance operator, which, up to negligible error terms, maps surface displace- 
ments to surface tractions of solutions of Lu = 0. This operator is introduced in 
Section [51 



6. Ray theory for systems 

The ray method solves a wave equation Au = by solving the equations F^j = 
for the amplitudes U-k of an asymptotic expansion of u. For j > 0, the 
equations F^j — are ordinary differential equations along rays, commonly called 
transport equations. 

The equation Fq = is the dispersion equation, Ao{x,d' (x))Uo{x) = 0. If Aq 
is scalar, then this becomes the eikonal equation Ao{x,9'{x)) = 0. Furthermore, 
F_i = with ([2^ is a transport equation for Uq which is solved by the method 
of characteristics. If Ao is not scalar, then it is not obvious how and F_i — 
lead to a useful transport equation for Uq. 

For systems of real principal type, introduced by iDenckeil (|l982[ ). there IS an 



efficient ray theory based o n the theory of Fourier integral operators. It applies 
to isotropic elastodynamics, Hansen fc Rohri j ( 2004 ). and, as we will show, to the 



propagation of Rayleigh waves along curved surfaces of general elastic media. We 
adapt this ray method to our setting. 

Let A he a (square) system of pseudo-difTerential operators with principal symbol 
Ao{x, ^). The set or zeros of the determinant det Aq is called the characteristic set 
of A (or of An). W e assume that the system A is of real principal type in the sense 



of iDenckeij (|1982| ). This means that there exists a matrix- valued symbol B{x,^), 
homogeneous of degree zero in the f variable, and a scalar symbol a{x, ^) such that 
B{x,^)Aq{x,£,) = a{x,£,)I holds with / denoting the unit matrix. Moreover, a is 
real- valued and its set of zeros equals the characteristic set, on which it vanishes 
simply, i.e. d^a{x,£,) ^ holds whenever a{x, £,) = 0. We call a a Hamilton function 
of ^0 (and of the operator A). In the special case where the determinant of Aq 
vanishes simply on the characteristic set, i? is a scalar multiple of the cofactor 
matrix, and ab = det Aq with b a nowhere vanishing function. 

Next w e collect some properties which follow from the real principal type as- 



sumption (Dencker, 19821) . The simple vanishing of a implies that the complement 
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of the characteristic set is dense in phase space. Moreover, BAq = al also AqB — al 
holds. This is clear where a ^ 0, and, by continuity, this holds in general. The 
following identities between range spaces and null spaces (kernels) are important: 

(27) range(_B) = ker(Ao), range(Ao) = ker(_B), where a = 0. 

The inclusion of the ranges in the null spaces is clear. Differentiating BAq — al — 
AqB in an appropriate direction, with the directional derivative denoted by a prime, 
we obtain 

BA'o + B'Ao^a'I = A'„B + AoB\ a' ^ 0. 

Now, equality in (|27|) follows from rank considerations. 

Fix a real-valued solution 6{x) of the Hamilton- Jacobi equation a{x,d'{x)) — 0. 
The integral curves (a;(s),^(s)) of Hamilton's canonical equations, x — d^a{x,£,) 
and ^ = —dxa{x,£,), satisfy ^(s) = 9'{x{s)) for all parameters s if this is true for 
at least one s. These integral curves in phase space are called bicharacteristics. 
The Poisson bracket {a, 6} is the derivative of a function b{x^^) in bicharacteristic 
direction, 

^b{xis),as)) = {a,b}ixis).m)- 
as 

The projections of bicharacteristics to space-time are the space-time rays x{s) as- 
sociated with 9. The vector field V{x) = d^a{x,9'{x)) is called the ray field. We 
denote by a dot the derivative in ray direction, 

tlixis)) = (d/ ds)Uix{s)) = {V ■ d^U){xis)). 

Next, we turn to the solution of Fq = and F_i = 0. Again we suppress 
arguments x and ^, and we assume evaluation of expressions at ^ = 9'{x). If 
AqUq = holds, then the leading term in the asymptotic expansion of Bf = BAu 
equals 

BF^i = Uo + (div V/2)Uo + iK^^Uo, 
where div V = Wx ■ V is the divergence of the ray field, and A'^^^ denotes the 
subprincipal symbol of A' = BA. To see this, observe that the right-hand equals 
the right-hand side in (1^^ when A is replaced by A' . Observe that A'g = al, and, by 
(P^ . A'^^^ = BAsuh + ^{^i ^o} when restricted to the nullspace of Aq. Therefore, 

= Uo + (div V/2)Uo + ^{B, A„}Ua + tBA^M- 

To solve F_i = we first solve, using the following lemma, the equations BF^i = 
and AqUo = 0, simultaneously. 

Lemma 6.1. IfU solves the differential equation 

(28) if + {div V/2)U+^{B,Ao}U + BW = 0, 

then AqU — holds along a ray if this holds in at least one point of the ray. 
Proof Multiply ^ by from the left to get 

Aoil + (div V/2)AoU + ^Ao{B, Ao}U = 
along rays. The derivative of Ao{x,^) along bicharacteristics is given by 

2^Ao = 2{al, Ao} = {AqB, Ao} - {Ao,BAo} = Ao{B, Ao} - {B, AojAo- 
as 



SUBSONIC FREE SURFACE WAVES IN LINEAR ELASTICITY 



13 



We obtain a homogeneous linear ordinary differential equation for AqU, 

-^{AoU) + idiYV/2)AoU+^{B,Ao}AoU = 0. 
as 2 

The assertion of the lemma follows from the uniqueness of solutions to initial value 

problems of this equation. □ 

Using the lemma, find a solution Uq =/= oi the transport equation 

(29) ilo + (divl//2)C/o + ^{B,Ao}Uo + iBA^M = 

such that AqUq — 0. Thus Fq = and BF^i = 0. The null space of B equals 
the range of ^o- Therefore, there exists C/_i such that F^i — —AqU-i. Set 
u — e'"^(J7o + (it^)^^t^-i) to obtain Fq = and F_i =0. To proceede by recursion, 
assume that u has been found such that F^j = and U-j+i = if j < fc. Then 

BF^k - U-k+i + (div V/2)U-k+i + lA^^^^U-k+i + M^-fc, 
where W-^k is a sum of derivatives of U-j, j < fc — 1. Using the lemma, solve 

il + (div V/2)U + Ao}U + iBA^^^U + BF^u = 0, A^U = 0. 

Replace U-k+i by U-k+i+U to obtain BF^^ = 0. (Note that F_j = still holds for 
j < k.) Choose U-k as a solution of + ^oC^-fc = 0. Then F_fe ==0, completing 
the recursive step of the construction of an asymptotic solution of Au — 0. 
For later reference, we summarize the result just obtained. 

Proposition 6.1. Let A pseudo- differential system of real principal type with 
Hamilton function a, and B(x,^)Ao{x,^) = o(a;,^)/. Let 9{x) with 9'{x) be a 
solution of the eikonal equation a{x,9'{x)) = 0. Denote by V{x) = d^a{x,0'{x)) the 
ray field, and div V = Va; • V its divergence. Let Uo{x) be a solution of the transport 
equation (j29p which satisfies 

Ao{x,0'{x))Uo{x)^O. 

Then there is a solution of Au = which is given by an asymptotic sum (1181) with 
leading amplitude Uo{x). 

7. Subsonic displacement boundary problems 

We return to elastodynamics. Given a subsonic phase 6 and a vector-valued 
amplitude W, each defined and smooth on the space-time boundary of the elastic 
body, we solve, asymptotically as — > oo, the displacement boundary problem 

Lu = in B, u = e^'^^'W at S. 

To construct u we use coordinates adapted to S as introduced in Section [51 Thus 
the boundary condition reads 

(30) uUa^o = e'^^'^'=°'^'''=''>W{x°,x\x^). 

The time variable is denoted t or x*', the dual variable is ^o- 

First, we define the subsonic region. Recall from ((TU)) the matrices Aj. Since the 
material parameters are allowed to vary smoothly with spatial position x, we have 
Aj — Aj{x, ri), and the acoustic tensor (1551) is given as follows: 

Lo{x,0 = (iAo{x) + ^3{Ai{x,rj) + Ai{x,rif) + A2{x,rj) - ^lp{x)I, 
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where 77 = (^^1,^2,0). We say that, at the surface point y = {x^,x'^), the covector 
C = {^o,£.i,£.2) is subsonic if Lq{x,£_)\x3=o is positive definite for real ^3. If ^ 
is subsonic, then so are — C and for s > 0. As in Section |31 we appeal to 
Proposition I A. 1 1 of the appendix, and we get, if C is subsonic at y, a factorization 

(31) Loix, = i^sl - QS(x, C))Ao{x){^3l - Qo{x, 0), 

where Qoix, C) depends smoothly on x and and the spectrum of Qo is contained 
in the upper complex half-plane. The factorization pip holds in a boundary layer 
< x^ < S. The matrices Qo{x, C) are uniquely determined. Using the uniqueness, 
it is seen that Qo satisfies Qo{x,sC) — sQo(x,(^) if s > 0. However, note that 
Qo{x, — C) 7^ —Qo{x, C) because of the spectral condition. It follows that Qo cannot 
be a polynomial in ^. 

The phase function 9 is assumed to be subsonic. By definition, this means that 
its derivatives C — 0' are non-zero and subsonic. Moreover, we assume that we have 
a zeroth order scalar symbol x(^)C)j supported in the subsonic region, such that, 
for 1^1 > 1, x(x,^) equals unity in an open conic neighbourhood Fi of the set F of 
{x,O,^ = se'ix),s>0. 

The pseudo-differential operator {xQo){x, D'), D' = (Do, fi, -D2), is defined, 
has the order one, and is tangential. A pseudo-differential operator is said to be 
tangential if it commutes with multiplication by a;^, or, equivalently, the operator 
pseudo-differentiates only with respect to the coordinates (a;°, x^, x^), and depends 
smoothly on as a parameter. 

The elastodynamic operator L is a second degree polynomial in D3 — —id^I 
with tangential differential operators as coefficients. Let Q be a tangential first 
order pseudo-differential operator, say Q = {xQa){x, D'), whose principal symbol 
is, on Fi, equal to Qo. By ([3T1) and the symbol calculus we have 

(32) L = iD3-P)AoiD3-Q)+R''D3 + R\ 

where P is the adjoint of Q, Pq — Qq, and is a tangential pseudo-differential 
operator of order at most j. Moreover, on Fi, the principal symbols R^ of R^ satisfy 

i?[J = Bo + id^Ao + AoQ-i + P-iAo - i ^ .(%Po)(a,Mo), 
Rl = B^- id^iAM - PoAoQ-i - P-iAoQo + ^ ^ X9^,^o)(a,MoQo)• 
Here we used (f23| and dH), and we have set 

(33) Bo{x) = [B'''% Bi(a;,ry) = 

Our next goal is to improve, by modifying Q and P, the factorization p2p so 
that the remainders Rj are negligible in a conic neighbourhood of F. To do so, we 
first eliminate P-i from the above equations for Rq, and we set i?Q = 0. We derive 
the following equation for Q_i: 

QqAoQ-i — AoQ-iQo 

^^^^ =BoQo + Bi-tAod3Qo + iY.idi,Q*o)Ao{d,,Qo)- 

This is a uniquely solvable Sylvester equation for AoQ-i because Qo and its Her- 
niitian adjoint have disjoint spectra. Having found Q-i, we solve for The 
symbols Q_i and P_i are only defined in Fi. Fix a new cutoff symbol xi which 
equals unity in a conic neighbourhood of F and is supported in the set where x — 1- 
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Replace Q by Q + {xiQ-i){^: D'), similarly for P. Then ([32l) holds with remainders 
which have, where xi — 1, order at most j — 1. This construction of reducing the 
order of the remainders can be iterated. In fact, one recursively solves equations 
QqX — XQo = Y for X — AoQ-k- Using asymptotic summation, a factorization 
([32| is obtained such that the remainders are negligible in a neighbourhood of F. 

Remark 7.1. The result ([32]) is a factorization of the elastodynamic operator into 
a product of pseudo-differential operators. The factorization does not hold every- 
where, but only microlocally in (some chosen subregion of) the subsonic region. 
Such a factorization is not possible using only differential operators, 

Our main application of the factorization p2p is the construction of solutions to 
displacement boundary problems. 

Lemma 7.1. Let9{x^,x^,x'^) be a subsonic phase function, W {x'^ , , x'^) a smooth 
amplitude. There is a solution u, 

(35) U - ^^u9(x°.x\x^) («w) ~' U-j {X°,X\X^,UJX'^), 

which solves Lu — asymptotically as w —> oo in x^ >0, and satisfies the displace- 
ment boundary condition (|30p . Moreover, D^u = Qu, and u decays exponentially 
in ujx'^ . 

The proof consists of finding an asymptotic solution (1351) of {D^ — Q)u — 
in x"^ > 0. To satisfy the boundary condition (pO| we require J/oU^^o — ^ f-i^d 
U-j\x3^Q = if j > 0. By the asymptotic expansion lemma for pseudo-differential 
operators, 

i{D3 - Q)u - Loe"^^ y°° {iLo)-^F-Ax\x^,Lox^), 

where, if we set C — 6'{x'^, x^,x'^), y = {x^,x'^, 0), 

F^j{x) = dsU-jix) - iQo{y,OU^,{x) + G-,{x). 

When deriving these formulas, the symbols Q-j of Q are replaced by their Taylor 
series in the x^ variable, and it is used that a power of x^ give rise to a power of 1 /w 
with the same exponent. The remainder term G-j is an expression involving U-£ 
for £ < j, and Go = 0. Wc recursively solve the equations F^j = for U-j imposing, 
in addition, the above-mentioned boundary conditions at the surface. Since iQo 
has its spectrum in the left half-plane, the solutions U-j{x) decay exponentially as 
x^ ^ oo. The factorization remainder WD^ + R^ is negligible in a neighbourhood 
of r. Therefore, (D3 — Q)u = implies Lu = 0. 

The assertion of the lemma still holds if pop is replaced by the more general 
displacement boundary condition u\r^3^Q = w, 

u;(a;0 , a;2) ^ ^ {iuj)-iW-j (x", x\x^). 

In fact, the proof of the lemma is readily adapted to cover this situation too. 



Remark 7.2. The methods of this section originate from the theory of elliptic bound- 
ary problems. Wloka et al. ( 1995h systematically use spectral factorizations of ma- 
trix polynomials in their treatment of such problems. 
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8. Existence of subsonic free surface waves 
As in (|15p we define, in the subsonic region, the surface impedance tensor Zq: 

(36) iZoiy, = Aoiy)Qoiy, C) + A,{y, r,). 

Here y = {x'^,x'^) = (a;\a;2,0), -q = (Ci,6) = (6,6,0), and C = (^o,??)- Since Qq 
and Ai are homogeneous of degree 1 in the variables ^j, so is Zq. The subsonic 
region is conic. Therefore, there is a function Coo{y, f]) > 0, homogeneous of degree 
one in rj, such that ( is subsonic if and only if Coo{y, ?/) > \£,o\ holds. If we freeze a 
point y on the surface and a unit horizontal propagation direction 77, \r]\ = 1, then 
Zq is the surface impedance tensor of a homogeneous half-space: Zo{y,C) = ^(c), 
where c — —6, and Coo (y, ''7) is the limiting velocity. This shows that ProDOsition l4.1l 
and the proof of the Barnett-Lothe Theorem apply mutatis mutandis to ZQ{y, C). 

Subsonic free surface waves correspond to solutions of the secular equation 
det^o(2/, C) = 0- We assume from now on that the secular equation has a zero 
6 < 0, necessarily unique and simple, for every surface point y and every direc- 
tion T] 0. By the implicit function theorem, there is a smooth function cji(y, rj), 
homogeneous of degree one in 77, such that < cr(jj, rf) < Coo(y, ?]), and 

(37) detZo-fo/i, h{y,0^^o + CR{y,v), %,C)>0. 

We call h the Hamilton function of the free surface wave problem. The phase 
function 9 = (p{y) — t is subsonic and solves h{y,6') = if and only if the eikonal 
equation Cfi{y,ip'{y)) = 1 is satisfied. The ray field is = (1, {dca/ dri){y , ip' {y))) . 
So we take time t as the parameter along rays, and we regard rays as the spatial 
curves y{t) which satisfy 

Note that the null-space N{y) of Zo(j/, —CR{y, »/), ?]), "q — (fi'iy), is one-dimensional. 
Denote by r > the distance to the boundary surface S. 

Theorem 8.1 (Existence of subsonic free surface waves). Let ip be a solution of 
the eikonal equation 

(38) CR{y,^'iy)) = l. 

There is a matrix-valued function H{y), which can be calculated algebraically from 
the derivatives up to second order of the elasticities, the material density, and the 
surface S, such that the following holds. Let Wo{y) satisfy, for every ray y{t), the 
transport equation 

(39) ^^Wo{y{t)) + ^dWV(y{t))Wo{y{t)) + H{y{t))Wo{y{t)) = 

and Wo{y{t)) G N{y{t)) for some t. Then Wo(2/) G N{y) holds for all y. The free 
surface boundary problem ^ has a non-zero solution 

(40) u(i,?/,r;w) - e"^(^(^)~*)y°° {iuj)-W-Ay,rio). 

The amplitudes U-j{y,ruj) decay exponentially as ruo — > 00, and the leading surface 
amplitude equals Uo{y,0) — Wo{y). 
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Let Q be the pseudo-differential operator of Lemma 17.11 with the defining cutoff 
symbol x chosen equal to unity in a conic neighbourhood of the set S which consist 
of {y,s(), where ( = {—l,(p'{y)), s > 0. Consider 

(41) T^[Tt], Tt^C^'''T%, 

and regard T as a multiplication operator. We define the surface impedance oper- 
ator Z by 

(42) Z=^i{AoQ + Ai+iT). 

Notice that Z is a first order pseudo-differential operator on the space-time bound- 
ary surface. The principal symbol of Z is —i{AoxQo + Ai) which equals the surface 
impedance tensor Zq in a conic neighbourhood of S. Let B{y, Q be the cofactor 
matrix of Zo{y, C) divided by the factor b{y, () of (I57|) . By Cramer's rule, 

(43) B{y,OZo{y,C) = h{y,Ol- 

Hence Z is (microlocally near S) a real principal type operator with Hamilton 
function h. 

We can now apply Proposition 16 . 1 1 to obtain a solution 



w 



of Zw = 0. (Notice that the transport equations can be solved with amplitudes 
which do not depend explicitly on time t.) From ([29]) we obtain the transport 
equation (I39p . where the coefficient matrix is given by 

(44) H =^{B,Zo}+tBZsub 

after evaluating at r/ = f'iy) and — —cniv, v)- Here B is defined in and ^sub 
is the subprincipal symbol of the surface impedance operator. Using Lemma 16.11 
it follows that Wo{y) G N{y). Using Lemma [7.11 and the remark at the end of its 
proof, we obtain a solution (j40p of the following displacement boundary problem: 
Lu = in r = a;'^ > 0, and u ^ w at — 0. Moreover, we get that the U-j decay 
exponentially as cjr — > oo. 

Recall the formula ([7]) for the surface traction r — [cTj'^] of the field u. Since 
D3U ~ Qu holds, we find that, in view of the definition of the surface impedance 
operator, the surface traction vanishes: 

—IT = AqD^u + Aiu + iTu = iZw = a.t = 0. 



This proves the existence of a free subsonic surface wave. In the next section we 
explain how H(y) is evaluated. 

Remark 8.1. The theorem says, in particular, that cn{y,r])/\ri\ is the wave speed 
of a subsonic free surface wave at the surface point y in the horizontal direction 
rj = ip'{y). Notice that an ansatz (HH)) will solve ^ only if the phase function ip 
satisfies the eikonal equation ([38)) . 
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9. Evaluation of the transport equation 

Before we discuss the evaluation of the coefficient matrix H{y), let us remark 
how the transport equation p9l) can be reduced to a scalar equation. Assume given 
a smooth unit field W{y) which spans the one-dimensional null-space N{y). Then 
the leading amplitude equals Wq — ipW with a complex-valued function "0(2/) which 
satisfies, along rays y{t), the scalar transport equation 

(45) -^i;+hdivV)i; + {W\-^W + HW)'ilj^O. 

The angular brackets denote the inner product defined by the metric tensor, con- 
jugate linear in the first slot. If the surface S is not plane, then the inner product 
will not be constant along rays, in general. 

We assume that the elasticity tensor, the material density, and the boundary sur- 
face are known, including first and second order derivatives. These data determine 
the transport equation. We make no attempt to derive analytical solution formulas. 
Rather we explain how the transport equation can be set up ready for numerical 
computation. To do actual computations, one choses coordinates y = {x^,x^) on 
S. The well-known ray coordinates are one convenient choice. For computations 
at interior points, the unique extensions as adapted coordinates, x = (x^, x^, cc^), 
are used. The data determine the acoustic tensor or principal symbol of the elas- 
todynamic operator, Lq. At every point of S, the acoustic tensor is viewed as a 
matrix polynomial in the variable conormal to S sd y. A factorization algorithm for 
matrix polynomials has to produce the factorizations (j3ip . Because of homogeneity 
it suffices to evaluate Qo{y,£,o,v) only at unit vectors rj. To actually compute Qo 
one may have to determine eigenvalues, eigenvectors, and in general also Jordan- 
Keldysh chains, of the matrix polynomial. The result of the computation also gives, 
by the impedance tensor Zq. 

First and second order derivatives of Zq are also needed. To calculate these, 
we procede as in the proof of Proposition I4.1f [e|) . The surface impedance tensor 
satisfies the Ricatti equation 

(Zo(C) - (ry))A-i(Zo(C) + iA,{7^)) = ^2(7;) - ^oP^, 

where we suppressed dependency on x from the notation. Denote differentiation 
with respect to some chosen coordinate or parameter by a prime. Differentiate the 
Ricatti equation, and obtain a linear equation 

(46) Q*X -XQo = Y 

for the complex 3x3 matrix X = Zq. The right-hand side Y is evaluated as an 
algebraic expression in Zq, Aq, Ai, p, and in the derivatives A'q, A'l, and p' . Since 
Qq and Qq have disjoint spectra, the Sylvester equation (|46p is uniquely solvable. 
Knowing Zq we also know the derivative Qq of Qo — AQ^{iZQ — Ai). Differentiating 
(|46|). we obtain an equation for the derivative X': 

Q*qX' - X'Qo = Y'- {Q'qYX - XQ'q. 

We conclude that first and second order derivatives of Zq can be computed from 
those of the elasticities and the material density by purely algebraic computations. 

Next we indicate how evaluate the Poisson bracket term (jUJ. Recall that bB is 
the cofactor matrix of Zq. Therefore, a first order derivative {bB)' is an algebraic 
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expression in the elements of Zq and Zq. So we can evaluate {bB,Zo}. Observe 
from dSll) that 

(47) 6(y,?7) =%detZo(2/,eo,?7) s^t h{y,Co,v) = 0- 

We calculate 

{B, Zo} = b~^{bB, Zo} + B{b-^I, Zo} = b-^{bB, Zq} - b-^B{bI, Zq} 

= b-^{bB, Zo} + b-^{bl, B}Zo + b-^{h, b}I. 

Here we have taken the Poisson bracket of (|33|) with bl to obtain the last equality. 
Hence 

{B, Zo}W = b-^{bB, Zo}W + b-^{h, b}W if W{y) e N{y). 

It is well-known that {h, /} vanishes on the set of zeros of h, if / does. Therefore, 
in view of (j47|) . the scalar Poisson bracket {h,b} remains unchanged at the zero set 
of /i, if we replace b by d^^g det Zq. 

By (|2T|) and (|23|) . the subprincipal symbol of Z is given by 

Z_i = AoQ-i+iT. 

We know already how to evaluate the second order derivatives of Zq. The equation 
(|34l) for AoQ-i is a linear equation of the form (j46| . The terms on the right- 
hand side of (|34|) can be evaluated. This is also true for the normal derivative 
dsQo- Indeed, the factorization ((3T|). and the Ricatti equation for Zq hold also for 
small a;'^ > 0, so that differentiation with respect to is possible. Furthermore, 
the right-hand side of (|34|) contains the terms Bq and Bi. These are defined in 
(p3| and in the displayed formula after ([7]). Notice that the formulas for Bq, Bi, 
and for T, defined in (j4T|) . contain Christoffel symbols. The curvature of S enters 
into the transport equation only through these expressions. Summarizing, we have 
explained how to evaluate H{y)Wo{y) in (15^ if, as we can assume, Wo{y) belongs 
to the null-space N{y) of the surface impedance tensor. 

The transport equation involves the ray field V{y) which depends on the eikonal 
(p{y). The zeros of the secular equation are simple, so they can be numerically 
computed in an efficient way. Differentiating 

= det ZQ{y,-CR{y, 7]), Tj) 

by the chain rule, the derivatives of cr are found as algebraic expressions in the 
derivatives of Zq. The solution of the eikonal equation for (p can be reduced to the 
solution of ordinary differential equations by Hamilton- Jacobi theory. This also 
allows the computation of the ray field V, the divergence divF, and the rays y{t). 

Remark 9.1. In the special case of a homogeneous body which fills a half-space, 
the tensors Qo and Zq are independent of surface points. So the Poisson bracket 
vanishes, and so does the subprincipal symbol, because the Christoffel symbols are 
zero too. The scalar transport equation (1451) reduces to 

-^V +{^divV + i lm{W\-^W))^j ^ 0. 

iXL Zi CI t 

If W is constant, this is the differential equation for the spreading factor. The 
imaginary term corresponds to the fact that a choice of W{y) is unique only up to 
a phase factor e*", a{y) real. In general, for an inhomogeneous, anisotropic body 
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with curved boundary, we have seen how to evaluate the coefficient {W\HW) of the 
scalar transport equation numerically. It is desirable, however, to also have a good 
understanding of the coefficients in (|45p . The real part, if positive, would lead to a 
damping factor. The imaginar y part gives rise to the B erry phase first observed by 



ry 

Babich in the early 1960's; see lBabich fc KiselevI (j2004l ) and the references therein 



Appendix A. Spectral factorization of positive definite matrix 

POLYNOMIALS 

The purpose of this appendix is to state and prove the spectral factorization 
theorem for self-adjoint matrix p olynomials which is fu ndamental to our approach 



to Rayleigh wave theory. Refer to iGohberg et al\ (|l982l ) for a comprehensive treat- 
ment of matrix factorizations. 

Let iJ be a finite-dimensional complex Hilbert space, dimH — n. Let A(s) be a 
quadratic polynomial with values in the space of linear operators on H, 

A{s) = Aos^ + (Ai + Al)s + A2. 

(We denote the adjoint by a star.) A number s e C is called an eigenvalue of the 
polynomial if A{s) is singular. The spectrum of A is the set (j{A) of its eigenvalues. 
Assume that, in addition, the polynomial is self-adjoint, i.e., A{s) = A{s)* holds 
for all s, and that Aq is positive definite. If s is an eigenvalue of A then so is s. If 
A has no real spectrum, then A{s) is positive definite for real s. 

In the following, we abuse language and often call linear operators matrices 
despite the fact that we do not fix a basis of H. 

The following is a special case of Theorem 11.2 in iGohberg et d\ (|l982h . 



Proposition A.l. Assume that A{s) has no real eigenvalues. Denote by a+ and 
(T_ the intersection of the spectrum of A{s) with the upper and the lower complex 
half-plane, respectively. There is a unique matrix Q with spectrum contained in the 
upper half-plane such that 

(48) A{s)^{s-Q*)A„{s~Q) 

holds for all s. If j+ is a closed Jordan contour which contains Oj^ in its interior 
and (T_ in its exterior, then 

(49) qI A{tY^At^ I tA{t)-^At 

Ji+ "'7+ 

holds. The integrals are non-singular. 

If the coefficient matrices Aj depend continuously or difi'erentiably on some pa- 
rameters, then, in view of (jlH]), so does Q. This follows from the integral formula 
(091). 

We need some preparations before we can give the proof of the propostion. 

Polynomial division of A{s) by s — Q gives A{s) = (s — Q-)Aq{s — Q) -\- F with 
Q-Aq -\- AqQ Ai -\- Al = and Q-AqQ -I- F = A2. A simple calculation shows 
that the remainder F is zero if and only if 

(50) AqQ^ + {Ai + ADQ + A2 ^ 

holds. The divisibility criterion (|50|) is known as the solvency equation. 
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The Stroh companion matrix of the polynomial A{s) is the linear operator N on 
which is given in block form by 



(51) 



N 



-A2 + AlA-'Ai 
A direct calculation proves the identity 



(52) 



sAq + Al 
-Ao 



-AlA-^ 



Ais) 
-sAa - Al 



where s = si. Obviously, the spectrum of A equals the spectrum of N. Moreover, 
{u,t)'^ is an eigenvector of N with eigenvalue s if and only if A{s)u = and 
t — {sAq + Ai)u hold. Set L ~ [l O] and R ~ [O /] . Passing to inverses in 
((52)) . we find that 

^-l _ r ^„ ^\^^-l ; 



(53) 



A{s)-' = L{s-N)-'R 



holds for all complex numbers s. 

Remark A.l. We call N the Stroh companion matrix of the pol ynomia l A(s) be- 
cause, in elasticity where n = 3, N equals Stroh's sextic matrix. IStrohl (|l962n . In 
that setting, and in the main part of the present paper, A(s) is the acoustic tensor 
at 77 + sv, where 1/ is a surface normal and rj a horizontal propagation direction. 
The Stroh matrix is a convenient example of a linearization of A{s), which means 
that ([53|) holds. The standard companion matrix of A{s) could however also be 
used for this purpose. 

The holomorphic functional calculus of N assigns to a function /, which is 
holomorphic in a neighbourhood of the spectrum a of iV, the matrix f{N) = 
§ f{s){s — N)~^ ds/27ri. The contour of integration must be chosen such that its 
winding numbers around the eigenvalues of iV in the support of / are equal to one. 
If ctq is a subset of the spectrum of N, and if /g = 1 and /o = in neighbourhoods 
of (To and a \ ctq, respectively, then Pq = fo{N) is a projector, called the Riesz 
projector associated with ag. 

By (|53l). if / is holomorphic in a neighbourhood of the spectrum of N, then 

(54) f{t)A{t)-^dt = Lf{N)R, 

where the contour of integration must be admissable for the functional calculus. 



Proof of Provosition \A.1\ Denote by P± the Ricsz projector of N associated with 
cr±. Clearly, P+P- = = P-P+, and P+ + P- is the unit matrix. The range Y± 
of P± is an A^- invariant subspace of H^, and the direct sum decomposition = 
Y+ © y_ holds. Denote by N± the restriction of A^ to y±. Clearly, ct(A^±) = a±. 

For sufficiently large p > 0, we denote by 7± the Jordan contour which consists 
of the segment [—p,p] and the semicircle pe^**, < i < tt. By ([5i|) we get the 
identity 

±LP±R = — i A(t)-^ dt = — f A(t)-^ dt. 

The last equality follows when letting p ^ 00. Since A{t)^^ is positive definite 
for t real, the integral defines a non-singular matrix. We infer that dim Y± = n. 
Moreover, LP+ has rank n, and there is a unique linear map : H ^ which 
satisfies L^LP+ = P+. Clearly, L" ~ P+L^ has rank n. 
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We prove that the proposition holds with Q = LNP+L^. Observe that = 
LN^P+LK Using dMl), we find that 

(55) — / s^A{s)-^ ds = LN^P+R = LP+R. 

Using Cauchy's integral theorem, we get 

(AoQ2 + (Ai + ADQ + A2)LP+R = 0. 

This implies the solvency condition (|50|) . Thus A{s) — {s — Q-)Aq{s — Q) holds 
with some matrix 

Observe that LP+L"^ equals unity. In view of the definitions, this implies that Q 
is similar to A^+. Thus the spectrum of Q lies in the upper half-plane. By ([5^ . we 
have det(Ao) det(s — N) ~ det(yl(s)), and therefore 

det(s - A^_) det(s - N+) = det(s - Q_) det(s - Q). 

This implies det(s — iV_) = det(s — Q-)- Hence the spectrum of Q_ equals the 
spectrum N^, which is contained in the lower half-plane. This means that the 
factorization obtained is spectral. 

If A{s) = {s — Q -)Ao{s — Q) holds with the spectra of Q and Q- contained in ct+ 
and (T_, respectively, then ([^^ follows by Cauchy's theorem. The integral on the 
left-hand side of (|49l) equals 2TTiLP+R which is non-singular. The uniqueness of the 
spectral factorization follows. Since a is self-adjoint, we also have the factorization 
A{s) = (s — Q*)Ao{s — Q*_). Since the spectra of Q and of Ql are contained in the 
upper half-plane, this proves that Q_ = Q* . □ 

Acknowledgement. The author thanks Aleksei Kiselev for valuable remarks and 
suggestions. 
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